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1. 


INTRODUCTION 


Boundary- layer  transition  is  a problem  of  significant  engineering  im- 
portance, which  at  the  present  time  is  not  well  understood.  For  design 
purposes  the  engineer  usually  must  use  empirical  correlations  of  experi- 
mental data  to  characterize  transition.  One  mechanism  which  leads  to 
transition  is  the  amplification  of  small  disturbances.  When  the  disturbances 
are  infinitesimally  small,  their  behavior  is  adequately  described  by  linear 
stability  theory,  although,  as  Gaster  has  shown,  there  are  some  ambiguities 
in  applying  solutions  of  the  Orr-Sommerfeld  equations,  which  assume  paral- 
lel flow  to  nonparallel  boundary- layer  flows.  Nonlinear  stability  theory 

extends  the  linear  theory  by  including  the  nonlinear  terms,  but  still  requires 

2 

the  disturbances  to  be  small.  (See  Stewartson  for  a recent  review  of  this 
subject.  ) All  nonlinear  stability  analyses  known  to  this  author  make  use  of 
the  parallel  flow  assumption  and  are  valid  only  in  the  vicinity  of  the  critical 
Reynolds  number. 

Two  earlier  numerical  studies  similar  to  the  present  work  are  those  of 
3 4 5 

DeSanto  and  Keller  and  of  Fasel  ’ . Each  of  these  is  a study  of  boundary- 
layer  stability  carried  out  by  integrating  the  complete  two-dimensional, 


Gaster,  M.,  "On  the  Effects  of  Boundary- Layer  Growth  on  the  Flow- 
Stability,  " J.  Fluid  Mech.,  Vol.  66,  1974,  pp.  465-480. 

> 

'Stewartson,  K.,  "Some  Aspects  of  Nonlinear  Stability  Theory,  11  Polish 
Academy  of  Sciences,  Vol.  7,  1975,  pp.  101-128. 

^DeSanto,  D,  F.  and  H.  B.  Keller,  "Numerical  Studies  of  Transition  from 
Laminar  to  Turbulent  Flow  over  a Flat  Plate,  n J.  Soc.  Indust.  Appl.  Math., 
Vol.  10,  1962,  pp.  569-595. 

^Fasel,  H.,  "Numerical  Solution  of  the  Unsteady  Na vie  r-Stokes  Equations  for 
the  Investigation  of  Laminar  Boundary  Layer  Stability,  " Proceedings  of  the 
Fourth  International  Conference  on  Numerical  Methods  in  Fluid  Dynamics, 


Springer- Verlag,  Berlin,  1974,  pp.  151-160. 

5Fasel,  H.,  "Investigation  of  the  Stability  of  Boundary  Layers  by  a Finite- 
Difference  Model  of  the  Navie r-Stokes  Equations,  " J.  Fluid  Mech.,  Vol.  78, 
1976,  pp.  355-383. 
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unsteady  Navier-Stokes  equations.  DeSanto  and  Keller  input  both  small  and 
large  amplitude  disturbances  into  the  boundary  layer  and  looked  at  the 
resultant  flow  downstream  of  the  disturbance  source.  Fasel  has  to  date 
concentrated  on  small  disturbances  and  thereby  investigated  the  effects  of 
the  nonparallel  boundary  layer  on  linear  stability.  (Fasel  has  indicated  in  a 
private  communication  that  he  has  also  carried  out  large  amplitude  calcula- 
tions which  are  qualitatively  similar  to  the  results  presented  herein.  ) Solu- 
tions to  the  three-dimensional,  unsteady  Navier-Stokes  equations  have  been 
carried  out  by  Orszag  ’ and  he  finds  a behavior  similar  to  transition. 

The  present  work  is  the  initial  phase  of  a study  whose  purpose  is  to 
develop  a computer  code  for  solving  the  Navier-Stokes  equations  (or  an 
appropriate  simplified  version)  in  the  Reynolds  number  range  of  transition 
and  to  study  the  various  stages  of  bounda ry - laye r transition.  The  nonlinear 
stability  of  the  Blasius  boundary  layer  is  investigated  by  simulating  the 
physical  situation  which  occurs  when  a vibrating  ribbon  introduces  distur- 
bances in  a boundary  layer.  The  solutions  have  been  obtained  by  solving 
both  the  two-dimensional  Navier-Stokes  equations  and  the  simpler  parabolized 
vorticity  equations.  The  detailed  formulation  of  these  equations  is  given  in 
Section  2;  the  spectral  method  used  to  solve  these  equations  numerically  is 
described  in  Section  3. 

2.  MATHEMATICAL  FORMULATION 

Consider  the  flow  of  an  incompressible  fluid  over  a semi-infinite  flat 
plate.  If  the  flow  is  assumed  to  be  two-dimensional,  then  the  unsteady 
Navier-Stokes  equations  may  be  written  in  terms  of  the  vorticity  u>  and  the 
stream  function  )/': 


Orszag,  S.A.,  "Turbulence  and  Transition:  A Progress  Report,"  Pro- 
ceedings  of  the  Fifth  International  Conference  on  Numerical  Methods  in 
Fluid  Dynamics,  Springe  r- Ve  rlag,  Berlin,  1976,  pp.  32-51. 
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+ *P  U>  - 4>  CJ  = v(u  +(J  ) 

t y x x y xx  yy 


(J  = Ip  + I p 

xx  yy 


(1) 

(2) 


The  independent  variables  are  x,  the  coordinate  parallel  to  the  plate;  y,  the 
coordinate  normal  to  the  plate;  and  time  t.  Equations  (1)  and  (2)  will  be 
solved  in  dimensionless  parabolic  coordinates,  defined  by  the  complex  trans- 
formation 


x + iy  = Xj 


< +i"(2R*1) 


1/2 


12 


(3) 


where  x.  is  a typical  distance  from  the  leading  edge,  and  R is  the  Reynolds 

1 . X1  . . 
number  based  on  the  freestream  velocity  and  x^.  The  dimensionless  time  T 

is  given  by 


T = tU  /x 
oo  1 


The  dimensionless  dependent  variables  are  defined  by 


(4) 


ip  - (Zi/u^xj)172^  = (zvu^yZg 


MS 

In  terms  of  these  variables,  the  Navier-Stokes  equations  are  as  follows 


(51 

(6) 


+ 1 


■ V* + w*  *Wizi 


(7) 


(8) 
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The  term  on  the  left-hand  side  of  Eqs.  (7)  and  (8)  containing  Rx^  is  negligible 


2 2 

compared  to  unity  for  rj  « 2£  R 


Since  £ = 0(1)  and  Rx  = 0(10' 


the 


,.  . . X1  1 
subject  term  is  important  only  for  very  large  values  of  T).  At  large  values 

of  rj , the  vorticity  is  exponentially  small.  Therefore,  the  term  in  brackets 
on  the  left-hand  side  of  Eqs.  (7)  and  (8)  may  be  set  to  unity  with  negligible 
effect  on  the  solution.  (A  comparison  of  numerical  solutions  with  the  sub- 
ject terms  retained  and  dropped  has  confirmed  that  this  term  is  unimportant.  ) 


Equations  (7)  and  (8)  are  fourth  order  in  both  £ and  r),  and  consequently 
require  two  boundary  conditions  on  each  boundary  of  the  solution  domain. 

Most  of  the  boundary  conditions  are  obvious  from  physical  considerations, 
but  the  downstream  boundary  conditions  are  difficult  and,  as  will  be  seen, 
may  complicate  the  numerics.  One  way  to  alleviate  this  problem  is  to  "para- 
bolize" Eq.  (7)  in  the  £ direction.  The  most  accurate  way  to  do  this  is  to 

substitute  Eq.  (8)  into  Eq.  (7)  and  drop  the  terms  of  0(R"2).  This  results  in 

^ 1 

a third-order  system  in  £,  and  Eq.  (7)  becomes 


^ = v +yy«  -gnw^\  + [<v>/(2«\) 


(9) 


As  stated  previously,  Eqs,  (7)  and  (8)  are  the  Navier-Stokes  equations; 

Eqs.  (8)  and  (9)  will  be  referred  to  herein  as  the  parabolized  vorticity  equa- 
tions. [Some  authors  neglect  the  entire  last  term  in  Eq.  (7)  and  term  the 
resulting  set  the  parabolized  vorticity  equations.  ] 


For  infinite  Reynolds  number  and  steady  flow,  the  two  equation  sets 
may  be  reduced  to  the  Blasius  equation: 


f + ff 

rjrj 


= o 


(10) 


In  many  cases,  the  solutions  of  interest  are  small  relative  to  the  Blasius 
solution.  Thus,  to  improve  numerical  accuracy,  the  inhomogeneous, 
nonlinear,  perturbation  equations  obtained  by  subtracting  out  the  Blasius 
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solution  are  actually  the  ones  solved.  However,  the  equations  are  presented 
in  the  usual  homogeneous  manner  to  keep  the  algebra  simple  and  not  obscure 
the  physics. 

The  equation  sets  will  be  solved  in  the  semi-infinite  region: 

0 < rj  < oo  (12) 

The  boundary  conditions  at  the  wall  (7]  = 0)  are  such  that  the  two  velocity 
components  are  zero: 


The  station  rj  = oo  is  upstream  of  the  body  in  Cartesian  coordinates,  and  the 
stream  function  and  all  its  derivatives  are  known  (except  for  exponentially 
small  terms).  The  condition  on  the  stream  function  is 

g/£  = f - ri-p  (14) 

whereyS  is  a constant  characteristic  of  the  displacement  thickness.  An 
additional  explicit  condition  on  g^  or  Q is  not  required,  as  the  condition  is 
imposed  by  the  form  of  the  numerical  scheme  which  is  discussed  in  Section  3. 

Upstream,  at  £ = 1,  the  two  velocity  components  are  specified.  For 
the  results  described  subsequently,  the  upstream  boundary  condition  con- 
sists of  a linear  combination  of  the  Blasius  solution  and  a time  periodic  solu- 
tion of  the  temporal  Or r-Somme rfe Id  equation 

(15) 


g = frj,  ■ + A Re  [0(/7)  exp(-i<j  r)] 

6 Blasius  r r 


g.  = fn i . - 2AaIm  [0(77)  exp(-itL>  T)1 

Blasius  ^ r ,J 
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where  0 is  the  Orr-Sommerfeld  solution  for  a real  wave  number  a,  and  U)  is 
the  frequency.  The  imaginary  part  of  0)  which  causes  0 to  grow  or  decay  in 
time  is  ignored  and  replaced  with  some  fixed  amplitude  A. 

The  Navie  r-Stokes  equations  require  two  downstream  boundary  con- 
ditions; the  parabolized  vorticity  equations  require  only  one.  Roache  in 
his  discussion  of  outflow  boundary  conditions  used  in  various  finite  difference 
schemes,  suggests  that  conditions  on  the  higher  derivatives  are  less  restric- 
tive. The  conditions  imposed  on  the  Navier-Stokes  equations  set  the  first 
and  third  normal  derivatives  of  g equal  to  zero  by  imposing 

= 0 (16) 

The  condition  imposed  on  the  parabolized  vorticity  equations  sets  the  third 
derivative  to  zero: 


= 0 


(17) 


The  ideal  downstream  boundary  conditions  have  no  upstream  influence; 
in  steady  flow  there  is  no  boundary  region  generated  at  this  boundary,  and  in 
unsteady  flow  there  is  no  upstream  reflection  of  the  wavelike  disturbances. 
Unfortunately,  the  conditions  (16)  and  (17)  are  not  ideal.  However,  because 
the  calculations  are  carried  out  at  a very  high  Reynolds  number,  the  down- 
stream viscous  boundary  region  that  occurs  when  the  Navier-Stokes  equations 
are  solved  is  a very  small  part  of  the  total  computation  domain,  having  a 
thickness  Reynolds  number  on  the  order  of  10.  By  the  same  token,  the 
reflected  wavelike  disturbances  which  occur  for  both  systems  of  equations 
have  a limited  but  larger  upstream  influence.  For  the  conditions  used  herein, 
it  was  found  that  the  reflected  wave  decayed  to  a negligible  amplitude  in 


Roache,  P.J.,  Computational  Fluid  Dynamics,  Hermosa  Publishers, 
Albuquerque,  New  Mexico,  1972,  pp.  154-161. 
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about  a wavelength.  (The  length  Reynolds  number  of  the  Orr-Sommerfeld 

4 

waves  used  in  this  study  is  about  2x10  . ) The  main  problem  with  the 
downstream  boundary  condition  effects  is  not  that  they  invalidate  part  of  the 
solution,  but  that  the  resulting  regions  must  be  resolved  or  the  calculation 
will  be  unstable.  Thus  the  boundary  region  can  set  the  resolution  required, 
and  indirectly  the  time  step,  and  can  thereby  substantially  increase  com- 
puter storage  and  running  time.  The  main  reason  for  using  the  parabolized 
vorticity  equations  is  that  the  very  thin  downstream  viscous  region  is  absent, 
and  thus  the  resolution  and  time  step  requirements  are  not  as  severe. 

Both  equation  sets  are  parabolic  in  time,  and  the  condition  imposed 
at  T = 0 is 


g Mf  = *f- 


Blasius 


(18) 


The  Blasius  solution  has  been  used  as  an  initial  condition  because  the  present 
study  is  concerned  with  perturbations  about  this  solution.  Other  initial  con- 
ditions have  been  used,  but  none  of  these  solutions  is  reported  here.  The 
Blasius  solution,  of  course,  is  not  an  exact  solution  to  either  Eqs.  (7)  and  (8) 
or  Eqs.  (8)  and  (9).  A solution  of  either  the  Navier-Stokes  equations  or  the 
parabolized  vorticity  equations  with  Eq.  (13)  as  the  initial  condition  and  no 
perturbations  [A  = 0 in  Eq.  (15)]  will  produce  a correction  to  the  Blasius 
solution  of  order  R , consistent  with  the  asymptotic  expansion  given  by 

g X J 

Goldstein  . These  corrections  to  the  Blasius  solutions  have  been  generated 

numerically  because  this  is  a convenient  problem  for  testing  numerical 

algorithms.  When  unsteady  perturbations  are  introduced,  they  are  much 
- 1 

larger  than  0(RX  ),  and  thus  the  resultant  solutions  are  unaffected  by  the 
higher  order  steady  part  of  the  solution. 


Goldstein,  S.,  Lectures  on  Fluid  Mechanics.  Interscience,  London,  1 4b0 
pp.  136-144. 
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3. 


NUMERICAL  METHOD 


; 


The  Navier-Stokes  equations  and  the  parabolized  vorticity  equations 
given  in  Section  2 have  been  solved  with  a "spectral  method"  (Orszag  and 
Israeli^):  the  dependent  variables  are  expanded  in  a finite  set  of  orthogonal 
functions  in  both  space  dimensions.  This  results  in  a coupled  set  of  ordinary 
differential  equations,  with  the  coefficients  of  the  orthogonal  functions  as  the 
dependent  variables.  These  equations  for  the  coefficients  have  been  numeri- 
cally integrated  in  time  with  various  finite  difference  methods. 

The  orthogonal  functions  are  not  chosen  arbitrarily,  but  rather  must 
satisfy  certain  requirements.  The  following  criteria  have  been  used  for 
selecting  the  orthogonal  functions  used  in  this  work:  (1)  the  function  set 

should  span  the  same  space  as  the  corresponding  independent  variables; 

(2)  since  derivatives  are  being  evaluated,  the  function  set  should  be  such  that 
derivatives  of  the  functions  can  easily  be  expressed  with  the  same  function 
set;  (3)  because  the  equations  are  nonlinear,  it  should  be  possible  to  map 
rapidly  the  coefficients  of  the  orthogonal  functions  into  a discrete  set  of 
physical  points,  and  vice  versa;  (4)  finally,  the  orthogonal  functions  should 
converge  rapidly.  For  the  present  problem,  a form  of  Chebyshev  poly- 
nomial is  used  in  both  dimensions.  For  example, 


,ri,T ) 


(19) 


where  T is  the  Chebyshev  polynomial  defined  on  the  interval  zero  to  one, 

and  n is  a scale  factor.  (The  properties  of  Chebyshev  polynomials  are 
r 


9 

Orszag,  S.A.,  and  M.  Israeli,  "Numerical  Simulation  of  Viscous  Incom- 
pressible Flows,"  Annual  Rev,  of  Fluid  Mech.,  Vol.  6,  1974.  pp.  281-318. 


t 
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or  Lanczos 


10 


1 1 


given  in  various  books;  see,  for  example,  Fox  and  Parker  or  Lanczos  ~.  I 
The  functions  used  in  Eq.  (19)  are  defined  on  the  proper  space,  and  there- 
fore meet  the  first  requirement.  The  derivatives  of  polynomials  or  exponen- 
tial polynomials  are  polynomials  of  the  same  type,  and  thus  the  second 
requirement  is  satisfied.  An  expansion  in  Chebyshev  polynomials,  such  as 
Eq.  (19),  is  a Fourier  cosine  expansion  on  a distorted  physical  space.  Thus 
the  mapping  back  and  forth  between  physical  space  and  Fourier  space  may 
be  carried  out  using  the  "fast  Fourier  transform.  11  (An  efficient  method  for 

performing  cosine  transforms  with  a standard  fast  Fourier  transform  routine 

12  13 

is  given  in  Cooley,  et  al.  . ) Orszag  has  shown  that  polynomials  such  as 
those  of  Legendre  and  Chebyshev  converge  much  more  rapidly  for  boundary 
value  problems  than  periodic  functions  such  as  sine  and  cosine  because  of 
the  Gibbs'  phenomena  at  the  boundary  associated  with  the  latter  class  of 
functions. 


All  dependent  variables  and  their  derivatives  may  be  expanded  in  a 
series  of  the  form  (19)  except  the  stream  function,  which  is  unbounded  at 
infinity;  in  this  case,  the  expansion  of  g - 77 ^ , which  is  finite  everywhere,  is 
used.  This  representation  of  g together  with  Eqs.  (8)  and  (14)  is  sufficient 
to  ensure  the  desired  exponential  decay  of  vorticity  at  infinity. 
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The  solutions  have  been  updated  in  time  with  both  explicit  and  implicit 
finite  difference  methods.  The  explicit  Euler  method 

a..(r+4r)  - a (r)  = 4tr  (r)  (20) 

has  been  used  as  well  as  the  following,  fully  implicit  method 

a..(T+  At)  - a..(T)  = 0.  54t[r..(t  + 4t)  + R..(T)]  (21) 

ij  ij  ij  ij 

where  a.,  are  the  coefficients  of  Q as  given  in  Eq.  (19),  and  R..  denotes  the 

ij  ij 

Fourier  transform  (or  coefficients)  of  the  right  side  of  either  Eq.  (7)  or 

Eq.  (9)  divided  by  2£.  The  evaluation  of  the  linear  terms  in  R..  from  the 

lJ 

known  expansions  of  g and  Q involves  some  simple  matrix  multiplications, 

the  details  of  which  may  be  obtained  from  Fox  and  Parker^.  The  nonlinear 

terms  are  obtained  by  Fourier-inverting  the  four  quantities  appearing  in  the 

two  nonlinear  terms,  carrying  out  the  indicated  multiplications  at  each 

discrete  point  in  physical  space,  and  Fourie r-t  rans  forming  the  sum  of  the 

nonlinear  terms.  This  is  an  efficient  way  of  evaluating  the  nonlinear  terms, 

but  it  does  introduce  some  aliasing  error  into  the  result.  Orszag  and 
. 9 

Israeli  call  this  a "pseudospectral  approximation"  (as  opposed  to  a spectral 
approximation,  in  which  there  is  no  aliasing  error).  They  suggest  that  the 
extra  computation  time  required  to  remove  the  aliasing  error  does  not  im- 
prove the  accuracy  enough  to  justify  the  added  computation  time;  test  calcu- 
lations performed  in  the  preliminary  stages  of  the  present  work  support  this 
conclusion. 

Since  Eq.  (21)  is  nonlinear,  it  must  be  solved  by  iteration  on  the 

quantity  R.^(t  + At).  As  a first  guess,  for  a..(r  + At)  the  Euler  difference 

equation  (20)  is  used.  Subsequent  corrections  to  a . . (r  + At)  are  obtained 

ij 

from 

a"j+1(T+  At)  - a.j(T)  = 0.  5At\  r".  (t  + At)  + R.dT)]  (22) 
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where  the  superscript  denotes  the  n^  level  of  iteration.  When  a..  agrees 
with  a . to  within  some  specified  tolerance,  the  iteration  has  converged. 

Thus,  the  difference  method  is  in  effect  a predictor-corrector  method,  with 
the  number  of  corrector  steps  determined  by  the  convergence  test.  Although 
the  iteration  on  the  implicit  difference  scheme  will  converge  and  is  stable  for 
relatively  large  time  steps,  the  time  step  which  minimizes  computer  running 
time  is  small  enough  so  that  convergence  is  obtained  in  one  to  two  itera- 
tions. For  the  Navier-Stokes  equations,  this  optimum  implicit  time  step 
is  a factor  of  four  or  five  greater  than  the  maximum  stable  explicit  time 
step;  for  the  parabolized  vorticity  equations,  the  optimum  implicit  step  is 
approximately  equal  to  the  maximum  explicit  step.  The  result  is  that  the 
implicit  method  is  about  twice  as  fast  as  the  explicit  method  for  the  Navier- 
Stokes  equations,  and  about  half  as  fast  for  the  parabolized  vorticity  equa- 
tions. To  minimize  computer  time,  the  implicit  method  has  been  used  to 
solve  the  Navier-Stokes  equations,  and  the  explicit  method  to  solve  the  para- 
bolized vorticity  equations. 

Equation  (20)  or  (22),  together  with  the  vorticity  boundary  conditions, 
completely  determines  the  vorticity  (or  equivalently  its  expansion)  at  the 
new  time  or  iterate  level.  However,  as  is  often  the  case  with  the  vorticity- 
stream  function  formulations,  'ome  of  the  vorticity  boundary  conditions  are 
not  known  explicitly.  In  the  present  case,  the  two  velocity  components  are 
specified  on  the  wall  and  the  upstream  boundary,  resulting  in  conditions  on  g 
and  its  derivative  normal  to  the  boundary. 

Two  different  methods  for  determining  the  vorticity  boundary  condi- 
tions have  been  used.  The  method  which  is  easiest  to  implement  and  which 
requires  the  least  computer  time  is  an  approximate  scheme.  With  this 
method,  the  additional  boundary  conditions  are  substituted  in  the  Fourier- 
transformed  version  of  Eq.  (8)  for  the  equations  governing  the  highest  fre- 
quency terms  in  that  dimension.  This  Poisson-like  equation,  containing  both 
the  Dirichlet  and  Neumann  conditions,  gives  the  new  value  of  the  Fourier 
expansion  of  g,  and  the  expansion  of  £2  that  satisfies  the  boundary  conditions 
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follows  from  Eq.  (8).  The  more  conventional,  second  method  has  greater 
accuracy,  but  also  requires  more  computation  time.  Having  solved  the  finite 
difference  equation,  one  simply  guesses  the  boundary  conditions  on  Q and 
solves  the  Poisson  equation  for  g with  Dirichlet  conditions.  In  general,  the  g 
so  obtained  will  not  satisfy  the  desired  Neumann  conditions.  The  error  in 
the  normal  derivative  of  g is  then  used  to  correct  the  boundary  value  of  Q , 
and  the  Poisson  equation  is  solved  again.  If  the  corrected  Q is  not  exact,  it 
may  be  necessary  to  iterate  this  process. 

The  advantage  of  the  first  method  of  treating  the  boundary  conditions  is 
that  all  boundary  conditions  may  be  imposed  simply  and  the  Poisson  equation 
need  be  solved  only  once  per  iteration.  The  disadvantages  of  the  approach 
are  that  numerical  instabilities  may  occur  and  also  that  the  boundary  condi- 
tion information  is  somewhat  artificially  concentrated  in  the  high  frequency 
part  of  the  expansion  spectrum.  By  experimenting  with  various  matrix 

. operators,  it  was  found  that  the  numerical  instability  problem  could  be  cured 

by  zeroing  a column  in  the  second  derivative  operator  (Fourier  space  opera- 
tor) for  each  extra  boundary  condition  in  that  dimension.  The  concentration 
of  boundary  condition  information  at  high  frequencies  is  controlled  if  the 
diffusion  is  large  enough.  In  the  rj  direction  the  diffusion  is  very  large,  as  is 
obvious  from  the  boundary-layer  nature  of  the  flow,  and  this  method  of  im- 
posing boundary  conditions  can  be  used.  The  situation  in  the  £ direction  is 
much  different;  diffusion  in  this  direction  is  only  important  on  length  scale 
characterized  by  a Reynolds  number  of  order  10.  Using  the  Chebyshev 
expansion  (which  acts  very  much  like  a highly  variable-mesh  finite -diffe  rence 

la 

scheme),  it  is  possible  to  resolve  lengths  of  this  order  near  the  two  bound- 
aries, where  boundary  regions  may  occur,  but  not  in  the  center  of  the  com- 
putational domain.  Thus,  in  the  £ direction  it  was  necessary  to  use  the 
second  method  of  imposing  the  boundary  conditions.  The  method  used  for 
obtaining  the  correct  value  of  the  boundary  vorticity  is  described  in  the 
following  paragraph. 


14  15 

Israeli  ’ has  considered  the  problem  of  finding  the  vorticity  bound- 
ary condition  at  a solid  wall  when  solving  the  vorticity-stream  function 
equations  with  a finite  difference  scheme.  He  proposed  a rapidly  convergent 
iterative  scheme  for  finding  the  vorticity.  This  iterative  scheme  and  its 
convergence  properties  are  based  on  the  fact  that  the  value  of  the  vorticity 
at  a point  on  the  boundary  is  primarily  a function  of  the  normal  velocity  at 
that  point  and  a weak  function  of  the  normal  velocity  elsewhere  on  the  bound- 
ary. Stated  another  way,  the  matrix  relation  between  the  vorticity  at  the 
grid  points  on  the  boundary  and  the  corresponding  normal  velocities  is  diag- 
onally dominant.  Even  though  the  boundary  condition  of  interest  here  is  at 
an  inflow  boundary  rather  than  a wall  and  the  operations  are  carried  out  in 
Fourier  space,  the  situation  is  much  the  same.  The  relation  between  the 
Fourier  coefficients  of  P and  g^  at  the  upstream  boundary  is  diagonally 
dominant,  and  a similar  iteration  will  converge  rapidly.  In  the  present 
problem  there  are  usually  17  modes  in  the  r]  direction  [n  = 16  in  Eq.  (19)]; 
therefore,  instead  of  iterating,  the  whole  matrix  is  stored  (the  matrix  is 
actually  only  15  x 15  because  of  the  reduction  associated  with  the  boundary 
conditions),  and  the  Poisson  equation  need  only  be  solved  twice  per  step. 


A brief  outline  of  the  technique  used  to  solve  the  Poisson  equation  will 
complete  the  description  of  the  numerical  method  used  to  update  the  solution 
in  time.  The  inversion  of  the  T]  operator  in  Eq.  (8)  is  accomplished  using 
the  tensor  product  method  of  Lynch,  et  al.  . This  requires  computing  and 
storing  the  eigenvalues,  the  eigenvectors,  and  the  inverse  of  the  eigenvector 

14 

Israeli,  M.,  "A  Fast  Implicit  Numerical  Method  for  Time  Dependent 
Viscous  Flows,  " Studies  in  Applied  Math,  Vol.  LIX,  1970,  pp.  327-34^. 

Israeli,  M.,  "On  the  Evaluation  of  Iteration  Parameters  for  the  Boundary 
Vorticity,  '*  Studies  in  Applied  Math,  Vol.  LI,  1972,  pp.  67-71. 

Lynch,  R.E.,  J.  R.  Rice,  and  D.  H.  Thomas,  "Direct  Solution  of  Pa  rtia  1 
Difference  Equations  by  Tensor  Product  Methods,  " Numerische  Mathematik, 
Vol.  6,  1974,  pp.  185-199. 
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array  associated  with  the  discrete,  Fourie  r- space,  second  derivative  opera- 
tor in  rj . Since  only  17  to  21  modes  are  used  in  this  dimension,  the  storage 
requirements  are  very  modest.  In  the  £ dimension,  it  is  of  interest  to  run 
rather  large  problems;  for  the  results  described  herein,  a maximum  of  129 
modes  was  used.  The  inverse  in  this  dimension  is  carried  out  by  observing 
that  the  second  integral  operator  in  the  £ dimension  is  zero,  with  the  excep- 
tion of  the  main  diagonal  and  the  second  diagonal  above  and  below  the  main 
diagonal  (Fox  and  Parker^).  Thus,  with  some  manipulation  of  Eq.  (8),  the 

£ inversion  may  be  carried  out  by  using  an  adaptation  of  the  tridiagonal 

7 

algorithm  described  in  various  numerical  analysis  texts,  e.g.,  Roache  . 
Therefore,  the  complete  inversion  of  the  Poisson  equation  in  the  largest  case 
involves  two  matrix  multiplications  of  arrays  sized  17  x 17  and  17  x 129,  and 
17  solutions  of  tridiagonal  problems  involving  129  unknowns. 

4.  NUMERICAL  RESULTS 

The  response  of  the  flat  plate  boundary  layer  to  a periodic  disturbance 
of  varying  amplitude  imposed  at  R^  = 10^  is  described  in  this  section.  The 

disturbance  is  a solution  to  the  temporal  Or  r-Somme  rfe  Id  equation  at  R - 

5 x 

10  , such  that  the  real  part  of  the  complex  frequency  is 

U)  = U>  X,  /U  = 13.  16  (23) 

r r 1 oo 

Most  of  the  computations  are  carried  out  in  the  Reynolds  number  range 
10  < R^  < 2.  5 x 10  . The  initial  condition  is  the  Blasius  solution,  and  at 
time  zero  the  Orr-Sommerfeld  disturbance  is  turned  on  at  the  upstream 
boundary  with  a frequency  U)^  and  an  amplitude  A.  This  results  in  a Tollmien- 
Schlichting  wave  being  propagated  downstream.  The  calculation  is  run  until 
the  solution  is  time  periodic  at  all  points  in  the  domain.  Experience  has 
shown  that  the  solution  at  a given  station  is  periodic  about  one  period  after 
the  disturbance  first  arrives  at  that  station. 
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Fasel**  has  published  calculations  similar  to  the  ones  described  herein 
in  which  small  amplitude  disturbances  are  introduced  into  the  boundary  layer. 
It  is  of  interest  to  compare  the  neutral  stability  Reynolds  number  and  the 
amplitude  behavior  as  a function  of  Reynolds  number  computed  by  Fasel's 
Navier-Stokes  code  with  the  present  parabolized  vorticity  results.  At  a 

constant  value  of  y corresponding  to  77  = 58.  Z/J  R in  the  present  coordinate 

x 5 

system,  Fasel  found  that  the  neutral  stability  Reynolds  number  is  1.  36  x 10  . 

This  is  in  good  agreement  with  the  minimum  of  the  curve  shown  in  Figure  1. 

At  small  amplitudes,  the  disturbance  is  linear  and  may  be  adequately  charac- 

4 

terized  by  a dimensionless  amplitude  ratio  as  is  done  by  Fasel  . To  make 
the  comparison  in  Figure  1,  Fasel's  dimensionless  amplitudes  have  been 
scaled  such  that  the  "neutral  stability  point"  falls  identically  on  the  curve 
shown.  The  other  points  are  in  good  agreement  with  the  curve. 

To  gain  confidence  in  the  validity  of  the  numerical  computations  when 
large  amplitude  disturbances  are  input  into  the  boundary  layer,  various 
self-consistency  checks  have  been  made.  Figure  2 shows  one  such  compari- 
son in  which  a Navier-Stokes  solution  is  compared  with  a parabolized  vorti- 
city solution.  The  agreement  in  both  phase  and  wave  shape  is  seen  to  be 
very  good.  Both  of  these  calculations  were  done  for  a Reynolds  number 
range  10  < R^  < 2. 5 x 10  , and  17  modes  in  the  t]  dimension  (n  - 16).  How- 

ever, the  time  step,  the  £ resolution,  and  the  time  difference  method  were 
different.  Even  though  twice  as  many  £ modes  were  used  in  the  case  of  the 
Navier-Stokes  equations,  experience  suggests  this  calculation  would  ulti- 
mately be  unstable.  Because  the  effective  cell  Reynolds  number  at  the 
downstream  end  of  the  computation  is  not  0(1),  the  viscosity  dominated 
waves  reflected  off  this  boundary  are  not  properly  dissipated.  Instead,  these 
waves  propagate  back  upstream  and  destroy  the  whole  calculation.  The 
solution  region  depicted  in  Figure  2 is  free  from  any  influence  of  reflected 
waves.  If  the  Navier-Stokes  resolution  is  increased  sufficiently,  the  re- 
flected wave  is  confined  to  a finite  region  near  the  downstream  boundary; 
the  parabolized  vorticity  equations  also  admit  a reflected  wave,  but,  since 
the  viscous  length  scale  is  absent,  the  requirement  for  very  high  resolution 
is  eliminated. 
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To  further  investigate  the  validity  of  the  solutions,  parabolized  vor- 
ticity  calculations  were  made  in  which  the  resolution  in  each  space  dimen- 
sion was  changed.  The  r]  resolution  was  studied  by  increasing  n from  16 
to  20.  The  agreement  of  the  two  parabolized  vorticity  calculations  was  better 
than  the  agreement  between  the  parabolized  vorticity  and  the  Navie r-Stokes 
calculations  shown  in  Figure  2.  The  resolution  in  the  £ dimension  was 

5 

studied  by  repeating  the  nominal  calculation  made  on  the  interval  10  < 

c 5 5 

R < 2,  5 x 10  by  shortening  the  interval  to  10  < R < 2 x 10  , In  both 

X x 

cases,  65  modes  were  used  to  represent  the  dependent  variables.  Again, 
the  agreement  was  comparable  to  Figure  2. 

The  numerical  results  presented  herein  consist  of  two  calculations, 

5 

both  made  with  the  parabolized  vorticity  code  on  the  interval  10  < R < 

5 X 

2.5x10  , with  65  modes  in  £ and  17  modes  in  T).  The  calculations  are 

identical  except  for  the  amplitude  of  the  Orr-Sommerfeld  solution  used  in  the 

upstream  boundary  condition.  In  the  small  amplitude  case,  the  velocity 

perturbations  are  everywhere  less  than  one-tenth  of  a percent.  The  large 

amplitude  disturbance  is  80  times  the  small  one,  and  nonlinear  effects  are 

clearly  present.  In  most  cases  it  will  be  convenient  to  compare  the  small 

and  large  amplitude  results  in  order  to  clearly  illustrate  the  nonlinear  effects. 

As  a point  of  reference,  the  Fourier  amplitude  of  the  small  amplitude 

input  disturbance  is  shown  as  a function  of  T)  in  Figure  3.  The  Blasius 

solution  is  also  given  to  show  the  relative  variation  of  the  curves.  Finally, 

the  Fourier  amplitude  of  the  disturbance  after  it  has  propagated  downstream 
5 

to  R = 2.  2 x 10  is  shown.  Figure  3 shows  that  the  peak  in  T]  space  is 
shifted  as  the  disturbance  is  propagated  downstream  and  amplified.  Subse- 
quent curves  will  show  the  variation  of  the  disturbance  as  a function  of  R 
at  77  " 0.  2,  which  is  well  down  in  the  shear  layer,  and  f]  - 1.0,  which  is  near 
the  peak  disturbance  point. 

Figure  4 shows  a comparison  of  the  large  and  small  amplitude  velocity 
perturbation  at  the  T]  location  inside  the  shear  layer,  as  a function  of  R . 

While  the  small  amplitude  disturbance  closely  resembles  a modulated 
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sinusoid,  the  large  amplitude  perturbation  is  significantly  distorted.  How- 
ever,  as  noted  previously,  the  disturbance  at  any  station  is  periodic  in  time 
with  the  period  of  the  upstream  disturbance.  This  fact  suggests  that  a Fourier 
series  in  time  will  allow  a clearer  interpretation  of  the  physics  than  curves 
such  as  those  in  Figure  4. 

The  Fourier  amplitude  of  the  perturbation  velocity  is  illustrated  in 
Figure  5 at  the  same  station  in  the  boundary  layer  as  Figure  4;  the  Fourier 
amplitude  at  a station  near  the  location  of  the  maximum  disturbance  is  shown 
in  Figure  6.  The  large  amplitude  curves  (A  = 0.08)  in  Figures  5 and  6 
show  that  up  to  = 1.  3 x 10~*  the  primary  mode  is  changed  only  slightly 
from  the  linear  behavior,  while  the  second  mode  grows  substantially.  Thus 
in  this  range,  the  nonlinear  effects  which  feed  energy  into  the  second  mode 
are  destabilizing.  Beyond  the  maximum  amplitude  point  of  the  second  mode, 
the  relative  stability  becomes  a matter  of  definition. 

The  perturbation  velocity  is  herein  defined  relative  to  the  Blasius 
and  therefore  the  mean  perturbation  is  not  zero  in  the  large  amplitude  case. 

The  mean  perturbation  shown  in  Figures  5 and  6 is  clearly  less  important 
than  the  second  harmonic.  The  higher  harmonics  are  even  smaller  than  the 
mean  and  are  therefore  not  shown. 

Figures  5 and  6 suggest  that  the  amplitude  of  the  second  harmonic 

relative  to  the  first  harmonic  is  a function  of  T].  Thus,  the  variation  of  the 

harmonic  components  with  rj  is  illustrated  at  two  values  of  the  Reynolds 

number  in  Figures  7 and  8.  At  = 1.  3 x 10  the  large  and  small  amplitude 

primary  mode  shapes  are  very  similar.  The  shape  of  the  secondary  is 

qualitatively  similar  to  that  of  the  primary,  but  the  peak  and  phase  reversal 

points  are  closer  to  the  wall.  The  mean  perturbation  velocity  is  positive 

near  the  wall  with  a smaller  negative  region  adjacent.  (Only  absolute  values 

are  shown  in  the  figure.)  Figure  8 shows  a similar  result  at  R = 2.  2 x 10” 

x 

At  this  Reynolds  number,  the  primary  disturbances  are  different  for  the  two 
values  of  A.  The  mean  contribution  again  has  a positive  region  near  the  wall, 
followed  by  a now  significant  negative  region. 
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It  is  also  of  interest  to  consider  the  phase  relationship  between  the 
primary  and  secondary  modes.  The  sine  of  the  phase  angle  for  both  modes 
is  plotted  versus  in  Figure  9.  The  phase  angle  has  been  designated 
£*x  + 0 only  to  call  attention  to  its  higher  order  behavior.  That  is,  if  Cl  (or  0) 
for  either  curve  is  a constant,  then  0 (or  Cl)  is  a weak  function  of  x.  The 
numerical  results  show  that  the  secondary  is  very  nearly  a spatial  harmonic 
of  the  primary.  This  agrees  with  the  usual  assumption  made  in  analytic 
studies  of  nonlinear  theory,  that  the  second  mode  is  a harmonic  of  the  pri- 
mary in  both  space  and  time.  It  is  interesting  to  note,  however,  that  (1)  the 
secondary  is  not  exactly  a spatial  harmonic  of  the  primary  mode,  (2)  the 
phase  relationship  varies  with  R^,  and  (3)  the  phase  relationship  appears  to 
correlate  with  the  growth-decay  behavior  of  the  second  mode  as  depicted  in 
Figures  5 and  6.  If  the  two  curves  were  spatial  harmonics,  then  the  negative 

peaks  of  the  secondary  would  coincide  with  the  positive  and  negative  peaks 

5 5 

of  the  primary.  Between  1.6  x 10  and  1.8  x 10  , there  is  clearly  a change 

in  the  relative  phase  of  the  two  modes;  this  change  is  well  correlated  with 

the  minimum  amplitude  point  of  the  secondary.  Although  this  point  needs 

further  study,  a tentative  finding  of  this  numerical  work  is  that  changes  in 

the  relative  phase  of  the  two  modes  are  correlated  with  the  energy  interchange 

and  hence  the  growth/decay  behavior  of  the  second  mode.  (Work  carried  out 

1 7 

subsequent  to  the  initial  submission  of  this  report  expands  upon  and  verifies 
this  conjecture.  ) 

5.  CONCLUDING  REMARKS 

There  are  four  main  results  or  conclusions  of  this  work.  First,  it  is 
possible  to  solve  the  two-dimensional,  unsteady  Navier-Stokes  equations 
(as  well  as  the  simpler  parabolized  vorticity  equations)  using  an  orthogonal 
function  expansion  (spectral  method)  in  both  space  dimensions  with  nontrivial 
boundary  conditions.  The  most  difficult  problems  arose  with  regard  to  the 
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Fig.  9.  Sine  of  the  Phase  Angle 
Secondary  Wave  at  77  = ( 


boundary  conditions.  The  first  problem  is  what  mathematical  boundary 
conditions  should  be  used  (particularly  on  the  downstream  side);  second, 
given  the  boundary  conditions,  how  should  these  conditions  be  implemented 
to  be  consistent  with  truncated  series  expansion  for  the  dependent  variables? 
The  results  described  herein  demonstrate  that  these  equations  can  be 
answered  in  such  a way  as  to  produce  useful  numerical  solutions.  However, 
further  work  is  currently  under  way,  and  it  appears  there  is  ample  room 
for  improvement  in  the  treatment  of  boundary  conditions. 

A second  important  result  of  this  work  is  that,  for  the  type  of  flows 
considered  to  date,  the  parabolized  vorticity  equations  provide  an  adequate 
model  of  the  flow,  in  that  these  solutions  are  in  agreement  with  solutions  to 
the  full  Navier-Stokes  equations.  This  result  is  important  because  it  has 
been  possible  to  generate  solutions  to  the  parabolized  vorticity  equations  in 
as  little  as  one-twentieth  the  time  required  to  solve  the  same  physical  prob- 
lem using  the  Navier-Stokes  equations. 

The  third  result  of  this  work  is  that  the  nonlinear  effects  can  be  de- 
stabilizing relative  to  the  linear  effects.  That  is,  the  boundary  layer  is 
more  unstable  to  large  disturbances  than  to  small  ones.  This  is  important 
since  previous  nonlinear  stability  analyses  apply  only  to  parallel  flow  and, 
in  addition,  are  series  expansions  about  the  critical  point.  The  present 
results  are  for  a nonparallel  boundary  layer  and  may  be  obtained  for  arbi- 
trary Reynolds  number  and  upstream  disturbance. 

The  final  point  is  that  the  Fourier-transformed  results  presented  herein 
represent  a beginning  which  should  lead  to  a better  understanding  of  the 
mechanisms  by  which  energy  is  nonlinearly  exchanced  between  various  wave 
modes.  The  nonmonotonic  behavior  of  the  second  mode  is  particularly 
interesting;  the  fact  that  the  growth/decay  behavior  of  this  mode  appears  to 
be  correlated  with  the  relative  phases  of  the  first  and  second  mode  is  also  a 
new  result.  Further  numerical  and  analytical  studies  will  be  required  to 
determine  the  full  significance  and  generality  of  these  final  results. 
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